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Abstract 

In this paper, we begin with the nonlinear Schrodinger/Gross-Pitaevskii equation (NLSE/GPE) for modeling Bose-Einstein con- 
densation (BEC) and nonlinear optics as well as other applications, and discuss their dynamical properties ranging from time 
reversible, time transverse invariant, mass and energy conservation, dispersion relation to soliton solutions. Then, we review and 
compare different numerical methods for solving the NLSE/GPE including finite difference time domain methods and time-splitting 
spectral method, and discuss different absorbing boundary conditions. In addition, these numerical methods are extended to the 
NLSE/GPE with damping terms and/or an angular momentum rotation term as well as coupled NLSEs/GPEs. Finally, applications 
to simulate a quantized vortex lattice dynamics in a rotating BEC are reported. 
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1. Introduction 

The nonlinear Schrodinger equation (NLSE) is a partial differential equation (PDE) that can be met in many 
different areas of physics and chemistry as well as engineering [1 2, 84 92, 158 i 1 1651 1 1 701 . The most well-known 
and important derivation of the NLSE is from the mean-field approximation of many-body problems in quantum 
physics and chemistry [ 158 1, especially for the study of Bose-Einstein condensation (BEC) lfTTlll581 . where it is also 
called the Gross-Pitaevskii equation (GPE) ||29lll 18||157I|158I . Another important application of the NLSE is for laser 
beam propagation in nonlinear and/or quantum optics and there it is also known as parabolic/paraxial approximation 
of the Helmholtz or time-independent Maxwell equations lfTll2l l84lll70l . Other important applications include long 
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range wave propagation in underwater acoustics Ill52l fl70l. plasma and particle physics [ 1701 . semiconductor industry 
11431 11461 . materials simulation based on the first principle 11011 . superfluids If5l l57ll . molecular systems in biology 

ma. 

In this paper, we consider the following time-dependent NLSE |29lll70l 

s 2 

isd t iff(t,x) = Vfy(f,x) + V(x)tff(t,x) + /(|i^(f,x)|V(f,x), xeR d , t > 0, (1.1) 

with initial data 

i/,( f = 0,x) = i/f (x), xeR d , (1.2) 

where i = V— I is the complex unit, t is the time variable, x e R d is the spatial variable with d — 1,2, 3, if/ := i//(t,x) 
is the complex- valued wave function or order parameter, V 2 = A is the usual Laplace operator and tAo := if/o(x) is a 
given complex-valued initial data. Moreover, < s < 1 is a dimensionless parameter. In most physics literatures, 
it is taken equal to s = 1; however in the semi-classical regime, we have < s <k 1, e being also called the 
"scaled" Planck constant. Function V := V(x) is a given real -valued external potential and its specific form depends 
on different applications and could also sometimes be time dependent lf27l l29l 11581 11701 . For example, in BEC, 
it is usually chosen as either a harmonic confining potential, i.e. V(x) = ^j- If27l |29l 11581 or an optical lattice 
potential, i.e. V(x) = A\ cosiL^x) + cos^y) + ^3 cos(Lj,z) in three dimensions (3D) with A\, Aj, A3, L\, L2 and 
L3 constants [27 291 112411125111581 . or a stochastic potential for producing speckle patterns II 1481 1 1581 : in nonlinear 
optics, it might be chosen as an attractive potential, i.e. V(x) = — H- HI [84]. The nonlinearity / := /(p) is a 
real-valued smooth function depending on the density p := \if/\ 2 e [0, 00) and its specific form depends on different 
applications [Q] |2l 127] |29j [T58l [1701 . In fact, when f(p) = A, with A a constant, then the NLSE ([□} collapses 
to the standard (linear) Schrodinger equation l92l 11651 . The most popular and important nonlinearity is the cubic 
nonlinearity ifTl 121 l29l fT58l fTTUll 

f(p)=fr, 0<p<™, (1.3) 

where f3 (positive for repulsive or defocusing interaction and negative for attractive or focusing interaction) is a given 
dimensionless constant describing the strength of interaction. Other nonlinearities used in nonlinear optics include 
the cubic-quintic nonlinearity f(p) - /3\p + fcp 2 fl~l [84] 11521 11701 and the saturation of the intensity nonlinearity 
/(p) = with p\i, j6i, Pi and Co given constants IfTl l84l 11521 11701 . In some applications, nonlocal nonlinearities 
can also be considered, and in this case f(p) = U * p which is a convolution, with U := U(x) the kernel |f29l l33l 
l44l l56l l65l 11851 . Specifically, for the Hartree potential case, U(x) = is the Green's function of the Laplace 
operator in 3D, and then, the NLSEdTTTJ can be re-written as the Schrodinger-Poisson system in 3D [56 38 1. Finally, 
we remark here that the NLSE <| 1 - 1 p is also called the Gross-Pitaevskii Equation (GPE) when the nonlinearity is 



chosen as the cubic nonlinearity as in ( 1.3 I, especially in BEC [27 291 11581 . Thus GPE is a special version of NLSE 
l271129i rrT51fT57lfT58l . 

There are many important dynamical properties of the solution i// to the NLSE (JTTTJ). Among them, here we 
mention several important ones that we will use to justify different numerical methods on whether they are still valid 
at the discrete level after the NLSE (JTTTJ is discretized by a numerical method. In fact, the NLSE (JTTTJ is a dispersive 
PDE and it is time reversible or symmetric, i.e. it is unchanged under the change of variable in time as t — > —t and 
taken conjugate in the equation. Another important property is time transverse or gauge invariant, i.e. if V — > V + a, 
with a a given real constant, then the solution \jj — > u/e~' a ' /s which immediately implies that the density p = \tf/\ 2 is 
unchanged. The NLSE (JTTTJ conserves many quantities. Among them, the mass (or wave energy in nonlinear optics) 
and energy (or Hamiltonian in nonlinear optics) conservation are given as ll29l [741 fT70l 

jV(f):=ll-A(f,-)ll 2 = f |<A(f,x)| 2 dx= f |<K0, x)! 2 ^ f \M*tdx := N(0), t>0, (1.4) 

J«J> jR d JW 

E(f):= f 
respectively, with 



-|V«A(r,x)| z + V(xMf,x)r +F(mt,x)\ z ) 



dx = E(Q), t>0, (1.5) 



F(p):= [ P f(s)ds, 0<p<+oo. (1.6) 
Jo 
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If there is no external potential in the NLSE ( 1.1 1, i.e. V(x) = 0, then the momentum and angular momentum are 



also conserved 129, 74l[l70l . In addition, the NLSE ( 1.1 1 admits the plane wave solution as ijr(t,X) = Ae' (k " a "\ where 
the time frequency «, amplitude A and spatial wave number k satisfy the following dispersion relation ll29l l74l 11261 
fTTOl 

W =^ + I/(|A| 2 ). (1.7) 
2 e 



In this case, if in ID with s = 1 and the nonlinearity is chosen as the focusing cubic nonlinearity ( 1.3 1 with p < 0, it 
also admits the well-known bright soli ton solution as ITO |2l |49l [84111521 1 1701 

iff B (t,x) = ^Lsech(A(x-vf-x ))e i( "^i ( " 2 ^ 2)f+9o) , x e E, t > 0, (1.8) 

where is the amplitude of the soliton with A a real constant, v is the velocity of the soliton, xq and 6q are the 
initial shifts in space and phase, respectively. Since the soliton solution is exponentially decaying for \x\ — > +oo, then 
the mass and energy are well defined and given by: N(i(/b) = —jf and EiipB) = 4$ + i^re- For other more dynamical 
properties of the NLSE ( |1.1[ ) such as dark (black and grey) solitons in ID under defocusing cubic nonlinearity, well- 
posedness and finite time blow-up, we refer to [29 741 11081 [1091 1 170 1 and references therein. 



For studying numerically the dynamics of the NLSE ( 1. 1 1, different numerical methods have been proposed in the 
literatures @[I2|22]|29]|33[23II3E^ l99TT7l9l [1201 [1301 fl44l [T56l [1731 [T82ll and references therein. 

The main aim of this paper is to review different numerical methods for solving the NLSE ( |1.1| ) numerically, compare 
them in terms of keeping different dynamical properties at the discretized level, stability and accuracy, discuss different 
absorbing boundary conditions for truncating the NLSE (JTTTJ onto a bounded computational domain, and extend them 
for solving NLSE with damping and/or angular momentum terms as well as coupled NLSEs. 

The paper is organized as follows. In Section[2j we review several popular numerical methods in the literatures for 
discretizing NLSE/GPE under simple boundary conditions and compare their advantages and disadvantages, discuss 
different absorbing boundary conditions for NLSE/GPE, and review some robust and efficient numerical methods for 
NLSE/GPE in the semiclassical regime. Extensions to NLSE/GPE with damping and/or angular rotating terms and 
coupled NLSEs/GPEs are presented in Sections [3] and |4] respectively. Numerical comparison of different numerical 
methods and some applications are reported in Section [5] Finally, we draw some conclusions and discuss future 
perspectives in Section|6] 

2. Numerical methods for NLSE/GPE 

2.1. Some popular numerical methods 

Here we present several popular numerical methods for discretizing the NLSE/GPE (JTTTJ and discuss/compare 
their properties including stability, accuracy, computational cost and how much properties of NLSE/GPE that these 
numerical methods can keep at the discretized level. For simplicity of notation, we shall introduce these methods in 
one space dimension (ID), i.e. d = 1 in \\.\\ . Generalizations to d > 1 are straightforward for tensor product grids 
and the results remain valid with modifications. For d = 1, the NLSE/GPE truncated on a bounded interval 
(a, b) with homogeneous Dirichlet boundary condition (|a| and b large enough so that the error due to the truncation 
is negligible) becomes 

e 2 

ied,t]/(t,x) = -—d xx i//(t,x) + V(x)(f/(t,x) + f(\i(/(t,x)\ 2 )t(r(t,x), a < x < b, t > 0, (2.1) 
i//(t, a) = (A(f, b) = 0, t > 0, (2.2) 
tf/d = 0, x) = tfr (x), a<x<b. (2.3) 

In this case, the mass and energy conservation collapse to the following 



N(t) := \W(t, 



r*b pb r*b 

I \ifr(t, x)\ 2 dx = I |<A(0, xfdx = I Wx)\ 2 dx := JV(0), t > 0, (2.4) 

Ja Ja Ja 



E(t) 



J a 



.2 



~-\dMUx)\ 2 + V(xWt,x)\ 2 + F(ty(t,x)\ 2 ) 



dx = £(0), t > 0. (2.5) 
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Choose a time step t := Af > and denote the different times as t n := m for n = 0, 1,2, . . .; choose the mesh size 
h :- Ax - ^j-, with 7 an even positive integer, and denote the grid points by xj :- a + jh, for j = 0, 1, ... ,J. Let iff* 
be the numerical approximation of i[f(t n , xj), for j = 0, 1, . . . , J and n = 0, 1,2, . . .; let iff" be the solution vector at time 
t = t n with components ((A")o<y<y- In addition, for any complex-valued vector <f> = (<f>j)o<j<j, we define the following 
finite difference operators 



t>j + i ~ <Pj 



< j<J-l; 



1 < j<J- 1. 



(2.6) 



Then the Cmnk-Nicolson finite difference (CNFD) method — in which one applies the second-order centered 
difference scheme for spatial discretization and the Crank-Nicolson scheme for time discretization — for discretizing 
( PI reads as 12911301 IgOl l79l [TT5lfTT9l[l96l 



t^ +1 - if,". 



£ 



where 



{D 2 x r +l ) j + {Dlr)] + [v(xj) +G(\r; l \ 2 , ir/)] 

F( Pl )-F(p 2 ) 



l// n . +l + l// n . 



1 <j<J- 1, n > 0, (2.7) 



G(p t ,p 2 ) = 



f(O Pl +(l-0)p 2 )dO, 0< Pl ,p 2 <oo. 



The boundary and initial conditions ( |2.2j >-(2.3 1 are discretized as 

0, 



n > 0; 



0<j<J. 



(2.8) 



The above CNFD method ( |2.7[ ) is time reversible or symmetric, i.e. it is unchanged if t — > -r and n + 1 <-> n; and 
its memory cost is 0( J d ) in cZ-dimensions with J unknowns in each direction. It conserves the mass (2.4 1 and energy 
d23) in the discretized level (29] [30] [31] EH [1151 , i e. 



/-i 



y-i 



;'=1 



7=1 



/V" := A l^'l 2 = h ^ l<l 2 = *2 ^ o( ^ )|2 := iV °' " " 

/-I 
7=0 



+ V{x j M" j \ 2 + F(\i^\ z ) 



i.n\2\ 



n > 0, 



(2.9) 



(2.10) 



which immediately implies that the CNFD method is unconditionally stable. In addition, it can be proven rigorously 
in mathematics that the CNFD method is second-order accurate in both time and space for any fixed s = sq — 0(1) 
Ol I3T1 l79l II 151 . However, it cannot preserve the time transverse invariant and dispersion relation properties 



of the NLSE/GPE (2.1 1 at the discretized level l29l[30l[3T| . Furthermore, it is an implicit scheme and its practical 
implementation is a little tedious. In fact, at each time step, one needs to solve a coupled fully nonlinear system which 
can be solved by a fixed point or a modified Newton-Raphson iterative method [7 , 8 9 29 30 3 1 1 and thus it might 
be very time consuming. In general, the computational cost per time step is much larger than 0(J d ), especially in 
2D and 3D. In fact, if the nonlinear system in (|2.7|i is not solved very accurately, e.g. up to machine precision, the 



numerical solution computed from (2.7 1 does not conserve the energy in (2.10 1 exactly |31 1. 

Due to the high computational cost of the CNFD method, in the literatures it is motivated for considering some 
semi-implicit methods in which a linear system is to be solved at every time step. Thus, the computational cost is 
significantly reduced, especially in 2D and 3D. One of these numerical methods is the relaxation finite difference 
(ReFD) method — in which one gives a name to the nonlinear term and makes a second-order approximation at time 
t„ — for the NLSE ((21) |6T1 



1 / n 

— Ik . 

?. V 7 



n+1/2 



+ It 



n-1/2 



)=/(l^| 2 ), 



i <;<■/-!, 



n > 0, 



u/". +1 - if/j 



"T 



(Dir +l ) j + (Dir)} + \ [v(xj) + iT" 2 ] « + , 



(2.11) 



where u 1 ^ 2 = iff" = ifroixj) f° r < j < J. The boundary and initial conditions (2.2H2.3I are discretized as 



(2.8 I. Similarly, this ReFD method is time reversible or symmetric, its memory cost is 0(J ) in li-dimensions with J 
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unknowns in each direction, and it is easier to be implemented than that of CNFD. It is an implicit scheme where, at 
each time step, only a linear system needs to be solved for example by the Thomas algorithm at the cost of O(J) in 
ID and by some iterative methods such as conjugate gradient (CG) or multigrid (MG) method lR)Tl[8"T l in 2D and 3D 
at the cost, in general, larger than 0( J d ). Thus it is computationally much cheaper than that of the CNFD method. It 
conserves the mass ( 1.4 1 at the discretized level as (2.9 1 |61 1, and thus it is unconditionally stable. In addition, it can 



be mathematically proven rigorously that the ReFD method is second-order accurate in both time and space for any 
fixed e - eo - 0(1) \61'\. Furthermore, for the NLSE with the cubic nonlinearity ( 1.3 1, this method also conserves the 
following discrete energy defined as 



E" :=hft. \ (Dtr) .f + h g [ v(XjWjl 2 + ^2 u n-m 

7=0 j=l 



n > 0. 



(2.12) 



Of course, it cannot preserve the time transverse invariant and dispersion relation properties of the NLSE/GPE (2.1 



in the discretized level I6TI . In addition, for general nonlinearity other than the cubic nonlinearity, no discrete energy 
conservation has been proven for the ReFD method yet. Another popular method is the following semi-implicit finite 
difference (SIFD) method — in which one uses the leap-frog scheme for the nonlinear term in time discretization — 
for the NLSE (JO) (SEEDED 



Zl LI 

2t 



n- 1 



£ 

"T 



1 < j < J - 1, n > 1; 



(2.13) 



or 



- yy n . 
2t 



n- 1 



£ 

~~4 



t// ,+1 + 1//" 



(D 2 x f +1 ) j + {^<r 1 )j\ + V(x } ) J 2 1 + fiWfW), 1 < j < J- 1, n > 1; (2.14) 



The boundary and initial conditions (2.2 1-( 2.3 1 are discretized as ( |2.8| l and the first step can be computed as [29, 30, 31 ] 

e 2 



IT 

S 



{D x f). + v( Xj )^ + mf)^ 



\<i<j-\. 



(2.15) 



Again, this SIFD method is time reversible or symmetric, i.e. it is unchanged if t — > -r and n + 1 <-> n — 1; its memory 
cost is 0( J d ) in ^-dimensions considering J unknowns in each direction. It is also much easier to be implemented 
than that of CNFD and ReFD, especially in 2D and 3D. It is an implicit scheme where, at each time step, only a 
linear system (whose coefficient matrix is time independent) needs to be solved, by the Thomas algorithm in O( J) 
operations in ID and by the direct fast Poisson solver via the discrete sine transform (DST) at the cost of 0( J d In J) 
in 2D and 3D. Thus, it is significantly cheaper than the CNFD and ReFD methods, especially in 2D/3D. In addition, 
it can be proven that the SIFD method is second-order accurate both in time and space for any fixed s = eo = 0(1) 
ll29l l30l [3T1 and a fourth-order (or six-order or even eighth-order) accurate method can be very easily constructed 
via the Richardson extrapolation technique 11811 . Of course, it is c onditi onally stable under the stability condition 

T < T„ 



5 max„< jSJ |V(x J )+/(| 1 A'jp)| 



(or r 



< 



£max <,<j |/(|iA"| 2 )| 



for (2.14i) for n > 0, which is independent of the 



mesh size h. It cannot preserve the time transverse invariant, dispersion relation and mass and energy conservation 



properties of the NLSE/GPE ^Jj in the discretized level 1291 1301 1311 . 

Another way to handle the nonlinearity in ( fTT) is to use the time-splitting technique ll29l [361 l39l l40l [451 IT241 1 1 25 1 
l66l 11201 11561 11731 1 1741 [182 1, i.e. from time t — t„ to time t = f n+I , the equation (2.1 1 is solved in two splitting steps. 
One solves first 



isd t if/(t, x) 



-dxxifrit, x), 



a < x < 



t > t„, 



(2.16) 



with the homogeneous Dirichlet boundary condition ( 2.2 1 for the time step of length t, followed by solving 



isd,tff(t, x) = V(x)(f/(t, x) + f(\tfr(t, x)| 2 )i/r(f, x), 



a < x <b, t > t„, 



(2.17) 
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for the same time step. Eq. (2.16 1 will be first discretized in space by the sine spectral method and then integrated 
(in phase space) in time exactly B29l [36l l39l l40l l53l . Multiplying (2.17 i by i/f(t,X) (conjugate of i//(t,x)) and then 
subtracting from its conjugate [29, 36, 39 40], we get for p(t, x) := |t//(f, x)\ 2 



d t p(t, x) — 0, t > t n , a < x < b, 



which immediately implies that the density p is invariant for any fixed x in the second splitting step (2.17 1, i.e. 



p(t,x) :— \ij/(t,x)\ 2 = \iff(t n ,x)\ 2 — p(t n ,x), t > t„, a < x < b. 



Plugging (2.19 » into (2.17l, Eq. (2.17 » collapses to alinearODE as 

ied t \f/(t, x) - V(x)y/(t, x) + f(\ijj{t n , x)| 2 )i//(f, x), a < x < b, t > t n , 
which can be integrated analytically as 

>f,(t, x) = e -<'[^)+/(W^)l 2 )]('->„)/ £ ^ t > ^ a < x < b 



(2.18) 



(2.19) 



(2.20) 



(2.21) 



In practical computation, from time t = t n to t = t n+i , one often combines the splitting steps via the standard Strang 
splitting 1 169 1 — which is usually referred to as time-splitting sine pseudospectral (TSSP) method | 

— as 

j 



= ^ e- hE ^ 2 (^\ sin( W (x,- - a)) = J] e^^' 2 ^), sin 



i=\ 

e -iT[V(Xj)+fW'f\ 2 )]/2e^(2) 



(2.22) 



1 <j<J- 1, n > 0, 



where ij/^ = if/"j + — for n > 0, i]r. — tf/o(xj) for < j < J, and (i/^ 1 ')/ for 1 < I < J — I, the discrete sine transform 



coefficients of the complex-valued vector ifr ' = {if/ Q , iff] 



, xf/ff with (/r^ 



nl 
b — a 



j-i 



7=1 



./ 



7=1 



Ijn 



J 



0, are defined by 



!</</-!. 



(2.23) 



Again, the above TSSP method is time reversible or symmetric, its memory cost is 0(J d ) in ^-dimensions with J 
unknowns along each direction. Its implementation is much easier than for CNFD, ReFD and SIFD. It is an explicit 
scheme since there is no need to solve a linear system and, at each time step, the computational cost is 0(J d In J). 
It conserves the mass (2.4 1 in the discretized level as (2.9 1 1 29j |36l [39] 00] |53), and it is therefore unconditionally 



stable. In addition, it can be rigorously proven that the TSSP method is second-order accurate in time and spectral 
order accurate in space for any fixed s = so — 0(1) [29, 621 |113|[139II151II1681I175I . In addition, it is time transverse 
invariant, i.e. when V(x) — > V(x) + a with a a constant, then i//" — > \p'\e~ lanT l e which implies that the density p" := |i/>"| 2 
is unchanged; it has the same dispersion relation as the NLSE/GPE (2.1 1, i.e. if V(x) = and i/fo(xf) = Ae lkx ' for 
< j < J, then the numerical solution from the TSSP method is (//" = Ae'^ kx '~ M " i with u, A and k satisfying the 
dispersion relation ( |1.7| i provided that J > k. However, it cannot preserve energy conservation properties of the 
NLSE/GPE (2.1 1 in the discretized level 11391 1301 . In fact, the numerical solution from the TSSP method actually 
coincides with the exact solution of a modified PDE at each time step. This shows the existence of a modified energy 
ll89l l97l l98l 11031 11051 preserved by the numerical scheme that is close to the exact energy if the numerical solution 
is smooth. However, some resonances may destroy the energy conservation on long time evolution (see Chapter 7 of 



[ 103|). Finally, for designing high-order, e.g. fourth-order accurate in time, time-splitting spectral methods, we refer 
to ll431l531fT74l[T75ll and references therein. 

In the literatures, in some applications where the solution of the NLSE/GPE is not smooth, e.g. with random 
potential ifTOl I%1 T1401 11481 . then the following time-splitting finite difference (TSFD) - in which the time-splitting 



6 



/ Computer Physics Communications 00 (2013) l- \23\ 



7 



is applied first and then the free Schrodinger equation (2.16 » is discretized by the CNFD method 
discretizing the NLSE/GPE EH M 02 EH) as 



_ e -iT[V(x 7 )+/(|^f | 2 )]/2 £ ^,(2) 



l <;'<•/- 1, 

0<j<J, n > 0, 



is used for 



(2.24) 



where = iffo(xj) for < j < J- Similarly, the above TSFD method is time reversible or symmetric, its memory 
cost is 0(J d ) in li-dimensions with J unknowns in each direction; it is easier to be implemented than CNFD, ReFD 
and SIFD. It is an implicit scheme where, at each time step, only a linear system (whose coefficient matrix is time 
independent) needs to be solved, which can be done in O(J) operations in ID and 0(J d In J) operations in 2D/3D. 
Therefore it is significantly cheaper than the CNFD and ReFD methods. The TSFD method is second-order accurate 
both in time and space for any fixed s — £q = 0(1) [29, 49|. In addition, it conserves the mass (2.4 1 at the discretized 
level as ( |2.9| > [29, 48 49 1 and thus it is unconditionally stable; it is also time transverse invariant. However, there is 
no energy conservation property of the NLSE/GPE ( |2.1| i at the discretized level ||29l |48l [49l . 

For comparison and convenience of the readers, we summarize the main physical and computational properties of 
the above popular numerical methods CNFD, ReFD, SIFD, TSSP and TSFD in Table[T] From this Table, we can see 
that the TSSP method shares the most properties as the original NLSE/GPE. In addition, it is very efficient and accurate 
as well as easy to be implemented in practical computations, especially in 2D and 3D [29, 36, 39, 40, 451 11741 . In fact, 
the TSSP method becomes more and more popular in practical computations, especially in the numerical simulation 
of Bose-Einstein condensation 1291 l36l l45l [1241 [1251 [1741 [1761 . 



Method 


TSSP 


CNFD 


SIFD 


ReFD 


TSFD 


Time Reversible 


Yes 


Yes 


Yes 


Yes 


Yes 


Time Transverse Invariant 


Yes 


No 


No 


No 


Yes 


Mass Conservation 


Yes 


Yes 


No 


Yes 


Yes 


Energy Conservation 


No 


Yes 


No 


YeQ 


No 


Dispersion Relation 


Yes 


No 


No 


No 


Yes 


Unconditional Stability 


Yes 


Yes 


No 


Yes 


Yes 


Explicit Scheme 


Yes 


No 


No 


No 


No 


Time Accuracy 


2 th or 4 th 


2 th 


2 th 


2 th 


2 th 


Spatial Accuracy 
Memory Cost 
Computational Cost 


spectral 
0(J d ) 
0(J d log J) 


2 th 
0(J d ) 
» 0(J d )f\ 


2 th 
0(J d ) 
0(J a 'log J)f\ 


2 th 
0(J d ) 
0(J d log /)[] 


2 th 
0(J d ) 
0(J d log J)f\ 


Resolution when < s <k 1 


9 




h = 0(e) 
T = O(s) 


h — o(e) 

T = O(s) 


h — o(s) 
t = o(e) 


h = o(s) 

T = o(s) 


h — o(s) 

T = O(S) 



Table 1. Physical and numerical properties of different popular numerical methods in the tZ-dimensional case with J unknowns in each direction. 



Remark 2. 1 If the homogeneous Dirichlet boundary condition ( |2.2| i is replaced by the periodic boundary condition, e.g. 
for BEC on a ring [ 29 39 , 40 , 49 1, the above numerical methods are still valid provided that we replace 1 < j<J—l 



by < j < J — 1 in all the numerical methods, and replace ip'i 



K 1 



by 0» 



uT/ 1 and ^ 

(2) _ 



the CNFD, ReFD and SIFD methods, the sine basis by the Fourier basis in the TSSP method, and ip^ 



(2) 



i/ff and 



(2) 

in the TSFD method, respectively. Similarly, if the homogeneous Dirichlet boundary 







4 Only for cubic nonlinearity 

'Depends on the solver for the nonlinear system 

6 lfd= 1,0(7) 

1 Ifd= 1,0(7) 

8 Ifd= 1,0(7) 

9 For cubic repulsive nonlinearity 

7 
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condition (2.2 1 is replaced by the homogeneous Neumann boundary condition, e.g. for the dynamics of dark solitons 
and their interaction [26 29 49 1, the above numerical methods are still valid provided that we replace 1 < j < J— 1 by 
< < J in all the numerical methods, and replace = = by = and tf/f + \ = t/r'J+J in the CNFD, 
ReFD and SIFD methods, the sine basis by the cosine basis in the TSSP method, and [j/^ = i//^ — by ifP^ = \f/^ 
and 4ffi\ - m tne TSFD method, respectively ll26l [29l [49l . Then the physical and numerical properties of these 
numerical methods are still valid. 



2.2. Other numerical methods 

In the literatures, there are many other numerical methods proposed for discretizing the NLSE/GPE. For example, 
Sanz-Serna 1991 1 1601 1 1631 1 164 1 proposed the following second-order finite difference (SSFD) method which is well- 
adapted to the computation of soliton-like solutions of NLSE/GPE 



V( Xj )+f 



2\ 



,n+\ 



1 <j<J- 1, n > 0, (2.25) 



and the boundary and initial conditions (2.2 1-( 2.3 1 are discretized as (2.8 1. Another one is the leap-frog finite difference 
(LPFD) method as 1150111901 



^7 ~ 



(D^.+[v(^) + /(i^r)]^ 



i <./■</ - 1, 



n > 1, 



(2.26) 



and the boundary and initial conditions (|2.2|i-(|2~3"j) are discretized as (|2.8[) and the first step can be computed as 



(2.15 I. In fact, in the physics literatures, for time discretization, the fourth-order Runge-Kutta (RK4) method was also 
used [4, 24, 70 69], which is not time symmetric. Thus, in general, it should be avoid to use RK4 for solving the 
NLSE/GPE. In some mathematics literatures, for space discretization, the 4th-order compact finite difference method 
has been used [138, 184] and the finite element (FE) method was also developed and analyzed (§]|3). However, since 
the computational domains for NLSE/GPE are usually simple and regular, the FE method for spatial discretization 
is usually not adapted in practical computations. Last but not the least, we want to mention that, in some physics 
papers, for solving NLSE/GPE in 2D or 3D, the alternating direction implicit (ADI) is first adapted to decouple the 
NLSE/GPE into dimension-by-dimension and then the SIFD or CNFD is used to discretize each NLSE/GPE in ID 
Eg] \UB EQl QUI E21 . For other more numerical methods for NLSE/GPE, we refer to Q [22] [24J EH [75] |76l 
E3l2ai82l2]|9l|95][lW3 and references therein. 

In general, these numerical methods have less good properties as those methods mentioned in the previous subsection 
and thus we omit the details here for brevity. 



2.3. Perfectly matched layers and/or absorbing boundary conditions 

In some situations for solving the NLSE/GPE (JTTTJ numerically, e.g. very long time dynamics and the potential 
is not a confinement potential in nonlinear optics ET1 [23), the simple boundary condition, such as homogeneous 
Dirichlet or Neumann or periodic boundary condition used in the previous subsections to truncate (or approximate) the 
original NLSE/GPE from the whole space problem to a bounded computational domain, might bring large truncation 
errors except that the bounded computational domain is chosen extremely large and/or time-dependent. Thus, in order 
to choose a smaller computational domain which might save memory and/or computational cost, perfectly matched 
layers (PMLs) (60) or high-order absorbing (or artificial) boundary conditions (ABCs) lTT2ll2Tll23ll59l ll02l|149l[T54l 
need to be designed and/or used at the artificial boundary so that one can truncate (or approximate) the original 
NLSE/GPE into a smaller bounded computational domain. Over the last 20 years, different PMLs II 1531 [T95 1 and/or 
ABCs lillBlQIlQIlQlDIiaiZni^ have been designed for solving the NLSE/GPE in the literatures. 

Here we simply review some of them for the completeness of this paper. 

In fact, PMLs were introduced by Berenger in 1994 for electromagnetic field computations |60| and they have 
been extended to NLSE/GPE by various authors recently 11153111951 . Again, here we present the idea in ID. Suppose 
that one is only interested in the solution of the NLSE/GPE \\.\\ with d — 1 over the physical domain (a, b). One 
introduces two layers with width Rq > at x — a and x = b and defines the following function 



S(x) :- I + Ro-(x), a :- a - R < x < b + R :- b, 
8 
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where R e C and cr(x) is a real-valued function which must be chosen properly and very carefully. The goal is to 
artificially damp the wave traveling inside the two layers to zero without modifying their dynamics in (a, b) by properly 
chosen R and cr(x) lfT2ll2Tl[T53lll95L Different choices have been proposed in the literatures |[T2lll53lll95l . Among 
them, the following choice works well for the NLSE/GPE (JTTTJ in ID with cubic focusing nonlinearity |[T2l|1531ll95l 



R = cr(x) = ^ 



(x — a) 2 , a < x < a, 

0, a < x < b, (2.27) 

(x-b) 2 , b<x<b, 



where 6 is a positive constant. Then, the NLSE/GPE (JTTTJ with d = 1 will be truncated (or approximated) as 
e 2 1 / 1 \ 

ied t \l/(t, x) = d x d x if/(t, x) + V(x)tff(t, x) + f(\ip(t, x)\ 2 )<p(t, x), a<x<b, t > 0, (2.28) 



2 S(x) \S(x) 

ifr(t, a) = iff(t, b) = 0, t > 0, (2.29) 
^(t = 0, x) = ip a {x), a<x<b. (2.30) 

The above problem can be discretized by CNFD, SIFD, ReFD and TSFD methods straightforward provided that we 
introduce the following finite difference operator to approximate gj-^d x (x^)5 A t/(x)) as 

\X=Xj 



1 -dx\-^TTd x d(x) 



S(x) \S(x) 



1 



2h 2 Sj 



1 - d j-i ~ \7r^— + c—) d J + 7T~ 



Sj-l/2 ' \ S J-l/2 S J+l/2l S j+l/2 



where S j - S(xj), S ;_i/2 = S(xj - h/2) and dj is the approximation of d(xj). This PML can be efficient and accurate 
for the NLSE/GPE in ID in some cases if one can choose 6 and Rq properly. Of course, extensions to 2D and 3D 
are a little bit more difficult and it might give bad numerical results or even not work in some situations when one 
chooses improper layer function S (x) and/or layer width Rq. Another way to design the PML is to introduce a damping 
potential which is also called as complex absorbing potential (CAP) or exterior complex scaling (ECS) 11341 11661 . 
i.e. choosing S(x) = 1 for a < x < b and replacing the real-valued potential V(x) by a complex-valued potential 



V(x) := V(x) - icr(x), where <r(x) is given in (2.27 i with 6 and <x two real constants to be determined based on the 
application. Then, the problem can be discretized by TSSP, CNFD, SIFD, ReFD and TSFD methods straightforwardly. 
Again, this CAP or ECS can be efficient and accurate for the NLSE/GPE in ID in some situations if one can choose 
6 and Rq suitably. The extensions to the 2D and 3D cases are direct. It might also lead to large inaccuracies in the 
calculations, more particularly for nonlinear problems. In general, PMLs are more robust from the accuracy point of 
view Ii l6ll . 

Another way is to find the Dirichlet-to-Neumann (DtN) map for the Schrodinger operator without and/or with 
external potential and nonlinearity at the boundary of the bounded computational domain by using the continuous 
and/or discrete Laplace transform |[]2[l3][ll[T5][l7][lJ3[T9l^ [155|. The DtN map is usually nonlocal and a 
time-dependent pseudo-differential operator, and thus its approximations which involve time fractional derivatives 
and integrals, are usually used in practical computations. These boundary conditions are also called as absorbing (or 
artificial or transparent) boundary conditions according to their properties. The goal is to avoid or at least to minimize 
the wave reflected back inside the domain while it should be outgoing. Since there are several very nice review papers 
on this part, we omit the details here for brevity and refer to 1 12 1 and references therein. Due to the nonlocal ABCs at 
the artificial boundary, one can usually solve the problem on the bounded computational domain by the second-order 
finite difference, such as ReFD or SIFD method. It is usually very hard to design spectral order method in space for 
the NLSE/GPE with nonlocal ABCs. 

For comparison of the performance and effectiveness of different PMLs and/or ABCs for NLSE/GPE and their 
applications, we refer to 1(17111351118611187 189| and references therein. 

2.4. Numerical methods in the semi-classical regime 

When < s <K 1 in the NLSE/GPE i.e. in the semiclassical limit regime, then the solution propagates 

waves with wave-length at O(e) in both space and time l39l |40l 1711 11141 [1241 11301 . The highly oscillatory nature 

9 
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brings severe numerical burdens in numerical computation of the solution of the NLSE/GPE ( |1.1[ ) 1 83l 1 1241 1 1 301 . In 
fact, in order to capture numerically "correct" physical observables such as density and current, different numerical 
methods request different mesh strategies (or e-scalability) for discretizing the NLSE/GPE (jTTTJ in the semiclassical 
limit regime, i.e. < s « 1. Based on the analysis and extensive numerical studies 1 39 . 40, 1241 11301 , it has been 
found that the e-scalability for TSSP is mesh size h = O(s), and time step r = (9(l)-independent of s and r = 0(e) 
for the linear Schrodinger equation and NLSE/GPE, respectively; for the CNFD, ReFD, SIFD and TSFD as well as 
many other finite difference methods, the mesh size is h — o(s) and the time step is r = o(e) ll39ll40l[T30lll44lll45l . 
Thus, among those numerical methods for discretizing the NLSE/GPE (JTTTJ directly, the TSSP method has the best 
resolution in the semiclassical limit regime. 

Of course, there are many other efficient and accurate numerical methods for solving the NLSE/GPE (JTTTJ based on 
the oscillatory structure of the solution. For the linear Schodinger equation, the Gaussian beam method [ 1041 1130111321 
I1331I137U 159I shows good resolution in spatial discretization, which in general requires h = O(yfs) and t = 0(1)- 
independent of s. For the details of the Gaussian beam method, we refer to a recent nice review paper [ 1 30 1 and 
references therein. However, in general, the Gaussian beam method cannot be extended to deal with the NLSE/GPE 
( |l.l| l due to the fact that the superposition is used in the method. We also mention here that there are some numerical 
methods in the literatures for solving the linear Schrodinger equation based on the Liouville equation through the 
Wigner transform, which can give good results in the linear case in the semiclassical limit regime ll83lll 141 1 1 301 . 

Another approach for dealing with the nonlinear (or linear) Schrodinger equation (JTTTJ in the semiclassical limit 
regime is through the Madelung (or WKB) expansion 1141117X1 . In the semiclassical limit regime, i.e. < s <SC 1, we 
denote the solution iff to the NLSE/GPE (JTTJJi as ifi E and use the following WKB ansatz (or Madelung transformation) 

GDH2BED 

iff e (t,x) = yjp e (t,x) exp^ ' 5 ^'^ j, xeR d , t > 0, (2.31) 

where the real-valued functions p e = |i^ £ | 2 and S e are the density and phase, respectively. Plugging the above WKB 
ansatz into the NLSE/GPE (jXTTJ and identifying real and imaginary parts, the NLSE/GPE is reformulated as a coupled 
system for density and quantum velocity v s = VS s , which is known as quantum hydrodynamic (QHD) system (made 
of a compressible, isentropic Euler system with Bohm potential) 1171117211 . This model was used in [22 , 72 , 83 901 11291 
to build an asymptotic preserving (AP) scheme which allows to compute numerical solution with time step and mesh 
size independent of s. However, the QHD system fails to represent the solution of the original NLSE/GPE near 
vacuum, i.e. p e — IT7T1I721 . and thus the scheme in [90 1 suffers from difficulties when the density p e vanishes in the 
domain and/or shocks or sharp changes happen in the QHD system. To overcome this drawback in the Madelung (or 
WKB) expansion for including vacuum, Grenier [ 1 16 1 introduced a modified Madelung (or WKB) expansion with the 
following ansatz 111 16117111 

iff e (t,x) = A e (t,x) exp( — (r ' (x> ), xeR d , r > 0, (2.32) 



where A e :- a e + ib E is a complex- valued function with a E and b E two real- valued functions. Plugging ( 2.32[ i into the 



NLSE/GPE <[LT} and collecting 0(1) and <9(e)-terms, one obtains a system for A e and S e (or v £ = VS E ) llTTllTTSl . 
Based on this formulation, an AP scheme was proposed in [73 1, which works up to the time when caustics happens. 
More recently, another AP scheme [63 129 1 was presented based on the Grenier's expansion with a regularized term, 
which can work after the caustics happens. For more details, we refer to 163, 73 j and references therein. 



3. Extension to NLSE/GPE with damping and angular rotation terms 

3.1. For damped NLSE/GPE 

As proven in l29ll74l[T70 'l. finite time blow-up may happen for the NLSE/GPE \1.1\ with a focusing nonlinearity, 



e.g. the cubic nonlinearity ( 1.3 1 with /3 < in 2D/3D. However, the physical quantities modeled by if/ := if/(t,x) do 



not become infinite, e.g. in BEC, which implies that the validity of (JXTTJ breaks down near the singularity. Additional 
physical mechanisms, which were initially small, become important near the singular point and prevent the formation 
of the singularity 11941 [108 109 1 . In BEC, the particle density p := \if/\ 2 becomes large close to the critical point 



10 



/ Computer Physics Communications 00(2013) 1 ]23\ 



11 



and inelastic collisions between particles which are negligible for small densities become important 11331 1941 11081 
109]. Therefore, a small damping (absorption) term is introduced into the NLSE/GPE (jTTTJ which describes inelastic 
processes [29 35 37l l94l[T08lll09l . We are interested in the cases where these damping mechanisms are important 
and, therefore, restrict ourselves to the case of focusing nonlinearity, i.e. f3 < 0, where [5 may also be time dependent. 
We consider the following damped NLSE/GPE ||29l [35l [37l IT081 1 09 1 



i d,if/(t,x) 



V 2 <A(f,x) + V(x) iff{t,n) + f(mx)\ 2 Mt,x) - i g(\^(t,x)\ 2 )^(t,x), 



t > 0, X € 



,K? = 0,x) = ^o(x), 



x e 



(3.1) 
(3.2) 



where g(p) > for p := |t/<i 2 > is a real-valued monotonically increasing function. 

The general form of ( |3.1) covers many damped NLSE/GPE arising in various different applications. In BEC, 
for example, when g(p) = 0, ( |3.1[ > reduces to the usual NLSE/GPE ( |1.1[ ); a linear damping term g(p) = 5 with 5 > 
describes inelastic collisions with the background gas; cubic damping g(p) = 5\p with 5\ > corresponds to two-body 
loss 1 162 1; and a quintic damping term of the form g(p) = 62P 2 with 62 > adds three-body loss to the NLSE/GPE 
( ]1 .1 j i 11621 . It is easy to see that the decay of the mass according to (3.1 1 due to damping is given by 



d r 

at J m 

Particularly, if g(p) = 5 with 5 > 0, the mass is given by 
N(t) 



m,x)\ 2 dx = -2 f g(mx)\ 2 M( t ,x)\ 2 dx<o, 



t > 0. 



= f \if/(t,x)\ 2 dx = e- 2dt N(0) = e- 26 ' f |iAo(x)| 2 dx, t > 0. 



(3.3) 



(3.4) 



Due to the appearance of the damping term, new ideas are needed to deal with them and different numerical 
methods have been presented in the literatures 1351 1371. In fact, the numerical methods such as TSSP, CNFD, ReFD, 
SIFD and TSFD methods for the NLSE/GPE (JTTTJ presented in the previous section can be easily extended to the 
damped NLSE/GPE ( |3.1[ ). For simplicity of notations, here we only present the TSSP method for ( |3.1| l with quintic 
damping term in ID, i.e. d = 1 and g(p) = 62P 2 with 62 > 0. From time t = t n to time t = t n+ \, the damped NLSE/GPE 
(3.1 1 is solved in two steps. One solves 

id,t//(t,x) = — d xx ip(t,x), a < x < b, t > f„, 



with homogeneous Dirichlet boundary condition for one time step of length t, followed by solving 

id,i{f(t,x) = V(x)tl/(t,x) + f(\i//(t,x)\ 2 )iff(t, x) - iS2\^{t,x)\ A ^(t,x), a < x < b, t > t„, 



(3.5) 



(3.6) 



for the same time step r. Again, Eq. ( |3.5| l is discretized in space by the sine-spectral method and integrated in time 

in for 

(3.7) 



exactly. For t e [f„, f n+ i], multiplying the ODE (3.6 1 by if/(t, x) and then subtracting from its conjugate, we obtain for 
p{t,x) := \Mt,x)\ 2 ||33H 

t > t n , a < x < b, 



d,p(t,x) - -2S2p i (t,x), 
which can be solved analytically as 

p{t n , x) 



pit, x) = 



V45 2 (f - t„) + p 2 (t n , x) 



t>t n , a < x < b. 



(3.8) 



Plugging (3.8 1 into (3.6 1, we get a linear ODE as 



idtif/(t,x) = 



V(x)+f 



p(t„, x) 



— i<>2 



P (t n ,x) 



^46 2 (t-t„)+p 2 (t n ,x)j 
which can be integrated exactly as for < s < r and a < x <b 

p 2 (t„,x) 



46 2 (t-t„)+p 2 (t„,x) 



il/(t,x), t > t„, a < x < b, (3.9) 



i[r(t n + s, x) — 4>{t n , x) exp 



V(x)s+ -p 2 {t n ,x)\n- 

4 46 2 s +p z (t n ,x) 

11 



p(t n ,x) 



^462U +p 2 (t„,x), 



\ \ 
du 



(3.10) 
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For cubic nonlinearity, the last term in the above equation can be integrated analytically. For a more general nonlin- 
earity, if it cannot be integrated analytically, one can use a numerical quadrature, e.g. the Simpson's rule, to evaluate 
it numerically 13511371 . Then, we can construct the second-order time-splitting sine pseudospectral (TSSP) method 
for the damped NLSE/GPE ( 3. 1 1 via the Strang splitting ||35ll37l ; the details are omitted here for brevity. 



3.2. NLSE/GPE with an angular rotation term 

In view of potential applications of BEC, the study of quantized vortices, which are related to superfiuid properties, 
is one of the key issues. Currently, one of the most popular ways to generate quantized vortices from BEC ground state 
is the following: impose a laser beam rotating with an angular velocity on the magnetic trap holding the atoms to create 
a harmonic anisotropic potential. Various experiments have confirmed the observation of quantized vortices in BEC 
under a rotational frame [ 3 , 5 , 6 , 29 70] 11421 . At temperatures T much smaller than the critical temperature T c , BEC 
in a rotational frame is well described by the macroscopic wave function iff :- \ff(t, x), whose evolution is governed by 
the following dimensionless GPE with an angular momentum term around the z-axis [5,6 291 1551150111071 1158 1: 

id,iff(t,x) = -ivV(f,x) + y(x)i/r(r,x)+y8|i/r(f,x)|V(f,x)-QL z i/r(f,x), xeR d , t > 0, (3.11) 



with initial data 



uj(t = 0, x) = iAo(x), x e R d , (3. 12) 



where d — 2 or 3 for 2D and 3D, respectively, f3 is a dimensionless constant describing the interaction strength, 
V := V(x) is a given real -valued potential which is usually chosen as a harmonic potential, Q is the dimensionless 
rotation velocity, and 

L z = -i(xd y -yd x ^J (3.13) 

is the z-component of the angular momentum operator L = (L x , L y , L Z ) T given by L = x A P, with the momentum 
P = -i V. The appearance of the angular momentum term means that we are using a reference frame where the trap is 
at rest. The above GPE is time reversible and time transverse invariant and it conserves the mass ( |1.4| i and the energy 
defined as 151 l6l l29l 1551 I5U1 [TUTl [B51 



E(t) 



JR' 1 



^|V<A(r,x)| 2 + V(x)|<A(f,x)| 2 + ^|<A(f,x)| 4 - Q<A(f,x) L#(f,x) 



dx = E(Q), t>0. (3.14) 



Due to the appearance of the angular momentum term, new difficulties are introduced in solving the GPE ( |3. 1 1 1 > 
for rotating BEC numerically. Several efficient and accurate numerical methods have been proposed in the literatures 
for discretizing it E21 EH [34J SI] |42l [50] [70l [ISS . In fact, the numerical methods such as TSSP, CNFD, ReFD, SIFD 
and TSFD methods for the NLSE/GPE (JTTTJ presented in the previous section can be easily extended to the GPE 
< |3 - 1 1 1 > with an angular momentum rotation. For conciseness, we only present here the TSSP method for ( |3.11| > in 2D, 
i.e. d — 2. 

From time t = t n to t — t n+ \, the GPE < |3 - 1 1 1 > is solved in two steps. One solves 

id t ip(t,x) = --Vtff(t,x) - QL,</<f,x), t > t„, (3.15) 

for one time step of length t, followed by solving the ODE 

i d t ifr(t, x) = V{x)(jj{t, x) + p\ij/(t, x)| x), t>t n , (3.16) 

for the same time step r. Similar to ( |2.17|i, Eq. ( |3.16| l can be solved analytically |29, 34 1. Some numerical methods 
have been presented for discretizing ( |3.15[ l. One numerical method is to adapt the polar coordinates (r, 6) in 2D 
such that the angular momentum rotation becomes constant coefficient, and then to discretize it in the ^-direction by 
the Fourier spectral method, in the r-direction by the second-order or fourth-order finite difference method or finite 
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element method, and in time by the Crank-Nicolson method. For more details, we refer to 



. Another method 



is to apply the Alternating Direction Implicit (ADI) method to decouple (3.15 1 into two sub-problems as 



id,iff(t,x) = ~d xx ilf(t,X) - iQyd x iJr(t,x), 



i d t tfr(t, x) 



d YV if/(t, x) + iQ.xd y (f/(t, x), 



t > t„, 
t > t„. 



(3.17) 
(3.18) 



Now the first problem (3.17 1 is constant coefficient with respect to x and the second problem (3.17 » is constant 
coefficient with respect to y, and thus they can be discretized in space by the Fourier spectral method and integrated 
in time exactly. Again, for more details, we refer to [29 50|. 

A different way to apply the time-splitting technique to the GPE ( |3.1 1[ ) is the following: from time t = t„ to time 
t = t n +u one solves 



1 Ixl 
id,if/(t,x) = --V 2 i/<f,x) + yf(r,x) - D.L z yj(t,x), 



x e 



for one time step of length t, followed by solving 



i d t i/f(t, x) = V(x) - ^- <Kr, x) + fWt, xWiftt, x), x e 



t > t n , 



t > t„, 



(3.19) 



(3.20) 



for the same time step r. Again, similar to (2.17 1, Eq. ( 3.20 1 can be solve analytically ||29ll4TII . Eq. (3.19 1 can be 
discretized in space by the generalized-Laguerre-Fourier spectral method and integrated in time exactly. One of the 
advantages of this method is that there is no need to truncate the original GPE ( |3.1 1[ > onto a bounded computational 
domain. Again, for more details, we refer to [29 [41]. 

Very recently, a simple and efficient numerical method has been proposed for discretizing the GPE (jXTTJ via a 
rotating Lagrangian coordinate 11201 l42l [1 121 . For any time t > 0, let A(t) be an orthogonal rotational matrix defined 
as 

cos(Qf) sin(Qf) 
- sin(£2f) cos(Qf) 



A(t) 



if d = 2, 



(3.21) 



and 



A(t) = 



( cos(Qf) sin(Qf) N 
- sin(Clf) cos(Qf) 
1 J 



if d = 3. 



(3.22) 



It is easy to verify that A l (t) = A T (t) for any t > and A(0) = /, with / the identity matrix. For any t > 0, we 
introduce the rotating Lagrangian coordinates x as 120111 121 



x = A" 1 (f)x = A r (f)x o x=A(0x, xeR d , 
and denote the wave function in the new coordinates as <p :- <p(t,x) 

<f>(t,x) := iff(t,x) = ^(*,A(0x), xeW 1 , t>0. 



(3.23) 



(3.24) 



In fact, here we refer the Cartesian coordinates (f , x) to as the Eulerian coordinates and (t , x) to as the rotating La- 
grangian coordinates for any fixed t > 0. Using the chain rule, we obtain the following cZ-dimensional GPE in the 
rotating Lagrangian coordinates without the angular momentum rotation term I20l l42l : 



id,<f>(t,x) = 



-\v 2 + W{x,t)+pW 



4>(t,x), 



x e 



t > 0, 



where W(x, t) = V(A(f)x) for x e W 1 and t > 0. The initial data (3.12 > can be transformed as 



4>(t = 0,x) = ifr(t = 0,x) = tff (x) := 4> (x) = 0o(x), 

13 



x = x e 



(3.25) 



(3.26) 
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Then the GPE (3.25 1 with the initial data (3.26 1 can be directly solved by the TSSP presented in Section 2. After 
obtaining the numerical solution (f>(t , x) on a bounded computational domain, if it is needed to recover the original 
wave function ijj(t, x) over a set of fixed grid points in the Cartesian coordinates x, one can use the standard Fourier/sine 
interpolation operators from the discrete numerical solution </>(f , x) to construct an interpolation continuous function 
over the bounded computational domain 1 64 , 1 67 1 , which can be used to compute if/(t, x) over a set of fixed grid points 
in the Cartesian coordinates x for any fixed time t > 0. For more details, we refer to [42 1. 



4. Extension to coupled NLSEs/GPEs 



In many applications, e.g. multi-components BEC Il25ll28l[29l [l58l and/or interaction of laser beams lfTll2ll55l 
11521 . coupled NLSEs/GPEs have been used for modeling different problems. For simplicity of notations, here we only 
consider coupled NLSEs/GPEs with two equations and cubic nonlinearity for two-components BEC and/or interaction 
of two laser beams. Extensions to coupled NLSEs/GPEs with more than two equations are straightforward. Consider 

e Ei Kama Ha nsa Hsu 



id,*//i(t,x) = 
id t uj 2 (t, x) = 



-V 2 + Vi(x) + J S 11 |*A 1 | 2 +ySi 2 |(A2| : 



with initial data 



y/i + A1//2, x e 
\f/2 + Atfri, x e 

yj i (f = 0, x) = i//f(x), ifr 2 (t = 0, x) = t/4 0) (x), X € I 



--V 2 + V 2 (X) + j e 21 |<A 1 | 2 +£22^21' 



t > 0, 
t > 0, 



(4.1) 



(4.2) 



Here (i/^, t/^) := (<Ai(x, f), ^(x, 0) is the dimensionless complex-valued macroscopic wave function, Vi(x) and V2CX) 
are two given dimensionless real-valued external potentials, fin, fin = Pi\ and P22 are given dimensionless real 
constants describing the interaction strength, and A is a given dimensionless real constant describing internal atomic 
Josephson junction in a two-components BEC (25] |27l |29l [158] [1791 [J80] [J92]| . This coupled NLSEs/GPEs conserves 
the total mass as |2j [27] EH [J58] [192) 



N(t) 

and the energy as 



= (|<Ai(f,x)| 2 + |^ 2 (*,x)| 2 )<*X = (\i/ff\x)\ 2 + \u>f\x)\ 2 )dx := JV(0), 



t >0, 



(4.3) 



E(t) := f [i (IV^I 2 + IV^I 2 ) + Vi(x)|«Ai| 2 + y 2 (x)^2| 2 + IfPiMi 



+ ^2im A +/3i2l<Ai| 2 |<A2l 2 + 2A ■ Re(^i^) 



dx = E(0), 



t >0, 



where Re(f) denotes the real part of the function /. In addition, if there is no internal Josephson junction in (4.1 
A — 0, the mass of each component is also conserved [25, 27 291 11581 [1921 



(4.4) 
, i.e. 



#!(*):= f \^i{Ux)\ 2 dx= f \ifrf\x)\ 2 dx, N 2 (t):= [ 



\Mt,x)\ 2 dx= f Wf(x)\ 2 dx, 



t > 0. (4.5) 



Different efficient and accurate numerical methods have been proposed in the literatures for discretizing the above 
coupled NLSEs/GPEs E3 W\ ED EU [1801 [Ml . In fact, the extension of the numerical methods TSSP, CNFD, 
ReFD, SIFD and TSFD for the NLSE/GPE ([□} presented in Section 2 is direct for the coupled NLSEs/GPEs gj}. 



For simplicity of notations, here we only present the TSSP method for (4.1 



From time t — t n to time t = ?„ + i, the coupled NLSEs/GPEs (4.1 1 is solved in two steps. One solves 

1. 



idtfrinx) = --V>i(f,x) +Aiff 2 (t,x), 



id t i/f 2 (t, x) 



L. 



-V 1 \p 2 (t,x) +Aifri(t,x), 
14 



t > f„, 

t > t n , 



(4.6) 
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for one time step of length t, followed by solving 



id,i/fi(t,x) = [Vi(x) + y 6n|iA 1 (?,x)| 2 + j3 l2 \^ 2 (t,x)\ 2 ] (Ai(f, x), t > t n , 
id t t/s 2 (t,x) = [v 2 (x) +/32i|iAi(f,x)| 2 + /? 22 |^ 2 (f,x)| 2 ]t/r 2 (f,x), f > f„, 



(4.7) 



for the same time step r. Similar to (2.16 1, Eq. (4.6 1 can be discretized in sp ace b y sine spectral method and then 
integrated in time exactly [25 , 28 , 29, 192|. Like for (2.17 1, for t e [t n , t„+{], Eq. (4.7 1 leaves p. 



|^| 2 andp 2 := |^ 2 | 2 



invariant [25 , 291 11921 . i.e. p,(f, x) = Pjfe, x) and p,(?,x) = p 2 (f„,x) for t„ < t < t n+ \ for any fixed x; thus (4.7 1 can be 
integrated in time exactly [25, 28, 29, 192 1. Then, we can construct the second-order TSSP method for the coupled 
NLSEs/GPEs (gl) v;'a the Strang splitting |25] El EH |54] [192) . We omit the details here for brevity. 

We remark here that the above TSSP, CNFD, SIFD methods have been extended to solve many other nonlinear 
dispersive partial differential equations arising from different applications. For details, we refer to the Schrodinger 
equation with wave operator [30, 32 1, the Schodinger-Poisson system Il39ll44lll93l , the Zakharov system 1147*1 l46l l79l 
II 15111311 . the Klein-Gordon-Schrodinger equations ll52l . and the Ginzburg-Laudan-Schrodinger equations [ 191 1 and 
references therein. 



5. Numerical comparison and applications 

5.7. Comparison of different numerical methods 

In order to compare the numerical performance and accuracy of different numerical methods, such as TSSP, 
CNFD, SIFD, ReFD and TSFD, for the NLSE/GPE ([LI}, we take d = 1, e = 1, V(x) = and f(p) = -p in fLl) , and 
the initial data iffQ in ( |1.2| l as 

iffoix) = A sech(A(x - x )) e i{vx+B °\ x e R, 



with A = 2, v = 1 and x Q = 9o = 0. Then the NLSE/GPE ( fTT) with ( fL2] l has the exact bright soliton solution ( [T~8|) , i.e 



ifr(t, x) - iffsit, x), with p = —\, A = 2, v=l and xo = #o = 0. In our computation, we take the bounded computational 

domain as the interval (a, b) with a — -15 and b — 20 and homogeneous Dirichlet boundary condition, which are large 

enough so that the truncation errors can be ignored. In order to quantify the numerical solution, we use the Z°°-norm 

of the error between the numerical solution df. and the exact solution dr(t„,Xj) as 

j J 

eUtn) := max \^(t n ,xj) - e'Z(t„) := max(|i/r(f„, Xj )\ - |^|), n > 0. (5.1) 

The functions e£> and e™ allow respectively to measure the phase error and the modulus error. 



h 


fco=0.5 


ho/2 


ho/4 


ho/8 


fco/16 


CNFD eZ, 

pin 


2.48 
1.89 


1.87E0 
6.98E-1 


4.28E-1 
1.46E-1 


1.03E-1 

3.52E-2 


2.57E-2 
8.68E-3 


ReFD el 

m 


2.48 
1.89 


1.87E0 
6.98E-1 


4.28E-1 
1.46E-1 


1.03E-1 
3.52E-2 


2.57E-2 
8.68E-3 


SIFD el 

pin 


2.48 
1.89 


1.87E0 
6.98E-1 


4.28E-1 
1.46E-1 


1.03E-1 
3.52E-2 


2.57E-2 
8.68E-3 


TSFD el 

pin 


2.48 
1.89 


1.87E0 
6.98E-1 


4.28E-1 
1.46E-1 


1.03E-1 
3.52E-2 


2.57E-2 
8.68E-3 


TSSP el 

m 


1.485 
1.408 


3.81E-4 
2.45E-4 


8.63E-9 
4.49E-9 


<lE-9 
<lE-9 


<lE-9 
<lE-9 



Table 2. Spatial error analysis on errors e£ m (/ = 5) of different numerical methods for the NLSE/GPE in ID under different mesh sizes h. 

To test the spatial discretization errors of the different numerical methods, we fix t = 10~ 5 such that the time 
discretization errors are negligible. Table [2] shows the spatial errors el m (t = 5) for different numerical methods under 
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different mesh sizes h. Similarly, to compare the temporal discretization errors of different numerical methods, we 
take h = 3.5 x 1CT 3 to get very small spatial discretization errors. Table[3]shows the temporal errors e%f(t = 5) for 
different numerical methods under different time steps t ll49ll . 



T 


T =0.1 


T /2 


T /4 


ro/8 


Tq/16 


CNFD 


e p 

1 CO 


2.62E-1 


6.65E-2 


1.64E-2 


3.88E-3 


7.28E-4 




m 
oo 


1.08E-2 


2.87E-3 


6.70E-4 


1.16E-4 


6.64E-5 


ReFD 


e p 


3.11E-1 


2.68E-2 


1.97E-2 


5.07E-3 


1.47E-3 




m 

C ' DO 


2.25E-1 


5.37E-2 


1.34E-2 


3.37E-3 


9.25E-4 


SIFD 


g" 


2.96E0 


6.46E-1 


1.91E-1 


6.62E-2 


2.61E-2 




/)7 

CO 


2.57E-1 


1.13E-1 


5.34E-2 


2.58E-2 


1.26E-2 


TSFD 


e p 

L - CO 


8.51E-1 


2.00E-1 


4.97E-2 


1.25E-2 


3.37E-3 




771 

OO 


4.10E-1 


9.03E-2 


2.20E-2 


5.49E-3 


1.44E-3 


TSSP 


e p 


5.17E-1 


1.40E-1 


3.57E-2 


8.98E-3 


2.25E-3 




m 
c co 


4.98E-2 


1.64E-2 


4.21E-3 


1.06E-3 


2.65E-4 



Table 3. Temporal error analysis on errors e^"'(f = 5) of different numerical methods for the NLSE/GPE jl.lj in ID under different time steps r. 

From Tabs. [2] & [3] as well as additional numerical results not shown here for brevity, it is clearly demonstrated 
that TSSP is spectral order accurate in space and second-order accurate in time while CNFD, ReFD, SIFD and TSFD 
are second-order accurate in both space and time. For more numerical comparisons, we refer to 11291 [49l |80ll and 
references therein. 



5.2. Applications 

In order to show numerical results for problems coming from applications, we consider the dynamics of the 
NLSE/GPE < |3 - 1 1 1 > with a rotation term in 2D starting from a quantized vortex lattice for a rotating BEC 11291 [3T1 |4T1 
|42][50], i.e. we take d — 2, ft = 1000 and Q. = 0.9 in ( |3.11| l. The initial datum in ( 3. 12 1 is chosen as a stationary vortex 
lattice which is computed numerically by using the method in [29, 51] with the above parameters and the harmonic 
potential V(x,y) = \(x 2 + y 2 ). Then, the dynamics of the vortex lattice is studied numerically by perturbing the 
harmonic potential from V(x,y) - \{x 2 + y 2 ) to V(x,y) = \{y\x 2 + yly 2 ) with (i) case I: y x = y y - 2, and (ii) case II: 
y x - 1-1 and y y = 0.9, respectively, at time t = 0. In the numerical simulation, we choose the bounded computational 
domain as [-32, 32] 2 with homogeneous Dirichlet boundary condition, mesh size h = 1/16 and time step t = 10" 4 . 
Figure [l] shows the contour plots of the density function p(t,x) :- |(A(f,x)| 2 displayed on [-17, 17] 2 at different times 
for cases I and II. 



6. Conclusion and perspectives 

Due to its massive applications in many different areas, the research on numerical methods and simulation as well 
as applications related to the nonlinear Schrodinger/Gross-Pitaevskii equations (NLSE/GPE) has been started several 
decades ago. Up to now, rich and extensive research results have been obtained in developing and analyzing efficient 
and accurate numerical methods for the NLSE/GPE and in applying them for simulating problems arising from many 
different areas, such as Bose-Einstein condensation (BEC), nonlinear optics, superfiuids, etc. Nowadays, numerical 
simulation has become a very important tool in theoretical and computational physics as well as computational and 
applied mathematics for solving problems related to NLSE/GPE and it has been used to predict and guide new exper- 
iments due to the advances in numerical methods and their analysis as well as powerful and/or parallel computers. Of 
course, due to its dispersive nature of the NLSE/GPE, when one is doing numerical simulation, he/she should choose 
the numerical method, mesh size and time step as well as the computational domain properly and carefully so that the 
numerical results reflect "correct" physical phenomena. 

The research in this area is still very active and highly demanded due to the latest experimental and/or technolog- 
ical advances in BEC, nonlinear optics, graphene, semiconductors, topological insulators, materials simulation and 
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t=0 t=1 t=4 




Figure 1, Contour plots of the density function p := \i//(t, x)| 2 for the dynamics of a quantized vortex lattice in a rotating BEC for case I (top two 
rows) and case II (bottom two rows). 

design, etc. It becomes more and more interdisciplinary involving theoretical, computational and experimental physi- 
cists and computational and applied mathematicians as well as computational scientists. Of course, there are still many 
important issues to be addressed. For example, extension and designing as well as analyzing new numerical methods 
for highly oscillatory nonlinear dispersive partial differential equations, especially coupled NLSE/GPE with other 
equations such as the Davey-Stewartson system |170|, Kadomstev-Petiashvili equations 11701 . coupled NLSE/GPE 

17 



/ Computer Physics Communications 00(2013) l- \23\ 



18 



with quantum Boltzmann equation for BEC at finite temperature [29 188], etc., are always welcome and highly de- 
manded. Another issue is to design efficient and accurate numerical methods and apply them for studying numerically 
NLSE/GPE with random potential flU ED EQl HH or stochastic NLSE/GPE JSH [SSI E21 ESI [J471 with applications 
and NLSE/GPE in higher dimensions, i.e. d > 3 for many-body problems in quantum chemistry and materials simu- 
lation and design. Here memory and computational cost might be extremely high and thus parallel computing and/or 
sparse grids as well as spatial and temporal adaptivity are very useful and essential. Last but not the least, efficient im- 
plementation of numerical methods and portable and readable programming are very important from the application 
point of view. Although there are plenty of research codes available in the community 167 68l[T50lll77l . a software 
package or tool box (with parallel implementation in 2D and 3D) is still very useful for applications. In summary, in 
order to solve challenging scientific and engineering problems and/or guiding and predicting new experiments related 
to NLSE/GPE, the interaction and close collaboration between computational and applied mathematicians and theo- 
retical and experimental physicists and chemists as well as computational and applied scientists becomes more and 
more essential in this research area. 
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